Correctness / robustness:
- GTF/GFF3: fold the `stop_codon` feature into the coding extent so cds_end
  includes the stop codon, matching Picard and the refFlat path. GENCODE emits
  stop_codon as a separate feature that CDS excludes; RefSeq/refFlat already
  include it.
- refFlat parser: enforce txStart<txEnd and exon/CDS within the transcript
  span, so the unchecked u32 subtractions in count_blocks/build_loci can't
  silently wrap on a malformed file — reject it instead.
- rna: saturating_mul on the `2 * end_bias_bases` short-transcript gate so a
  pathological `--end-bias-bases` can't wrap it and panic the coverage slices.
- gene-model format auto-detect: choose GTF vs GFF3 by which key/value
  separator appears first, so a GTF value containing '=' isn't misread as GFF3.
- gene-model parse errors now carry the file path and detected format.

Cleanup / placement:
- Move overlap_inclusive/encloses_inclusive out of gene_model (rna is the only
  caller) into rna.rs and drop them from the public API.
- Move the MateSpan helper struct out of the module-function section.
- intervals::read_file_contents back to private, and fix its stale doc: it is
  not shared with the gene-model reader (which streams via fgoxide).
- Wrap the rna `--help` long_about to house style; mark the unused FG_PINE
  brand colour `#[allow(dead_code)]`; fix the module doc that overstated what
  `--exclude-duplicates` drops.
- CHANGELOG: correct the rna output-file list (drop the non-existent
  .rna-coverage.txt; add the biotype/gene-expression charts).

Tests: rna-through-multi integration test; stop_codon extension (GTF + GFF3);
refFlat span validation; contig reconciliation; format auto-detect '='.

Deferred to follow-ups (#20, #21): reduce rna per-read CIGAR allocation;
evaluate the noodles gtf/gff parsers. Skipped as theoretical: the split-MC
insert-size dedup edge case.

GFF3 percent-encodes reserved characters (; = & , % tab newline) inside
attribute values; our hand-rolled parser left them encoded, so a gene name or
Name containing an encoded comma/semicolon came through mangled. Decode values
at extraction (gene_name / Name / gene_id / transcript_id / ID, biotype, Parent,
and the RefSeq region accession->name mapping). Parent is split on the raw ','
separator first, then decoded, so an encoded comma (%2C) inside one id is not
treated as a separator.

percent-encoding was already in the tree transitively; promoted to a direct dep.

This closes out the noodles-gtf/gff parser evaluation. A full switch to noodles
was ~3.5x slower on the gene-model parse (8.4s vs 2.4s on GENCODE v29 GTF, since
noodles builds a full percent-decoded attribute map per line) and left all the
source-specific grouping / biotype-key / RefSeq-accession logic ours -- so
percent-decoding, the one real gap noodles would have closed, is added directly
to the hand-rolled parser instead.

bases_overlapping_exons was the single largest self-time item in the
--threads 3 hot-path profile (~6% self). It summed the overlap of every
alignment block against every exon of a transcript — an O(blocks×exons)
product — even though a read's blocks and a transcript's exons are both
sorted by start and internally disjoint, so each block can only overlap a
short contiguous run of exons.

Replace the nested product with the same partition_point + break
merge-walk already used by overlap_with_intervals: for each block,
binary-search the first exon that can overlap and scan only until an exon
starts past the block. O(blocks + exons) with early termination.

Output is byte-identical on the 80M-read ENCODE BAM (all four rna TSVs).
Instructions retired drop ~1.5% (~600.4B -> ~591.2B, consistent across
runs); wall time on this laptop is thermal-noise-dominated and can't
resolve a change this size.

classify_bases and tally_insert_size each called rec_alignment_blocks on
the same record independently, so every paired read walked and allocated
its CIGAR alignment blocks twice. Parse them once at the top of accept()
into a local Vec and pass &[(u32,u32)] to both; classify_bases no longer
needs the record at all.

The block Vec is deliberately a local passed by reference, not a field on
the collector: a self-owned scratch buffer read in a loop that also
mutates other self fields forces aliasing-conservative reloads (measured
~3.5% regression previously), whereas a provably non-aliasing local does
not.

Output byte-identical on the 80M-read ENCODE BAM. Instructions retired
drop ~2.4% (~590.2B -> ~575.8B, consistent across runs).

The MC (mate CIGAR) tag was parsed up to three times per read: once in
classify_ribosomal and once in tally_strand (both via mate_alignment_end,
each of which copied the tag bytes and re-walked the CIGAR only to take
its reference end) and a third time in resolve_mate_span for insert size.
For a paired RNA read that path meant ~3 tag copies + 3 CIGAR parses +
their allocations, every read.

Parse it once at the top of accept() into a MateAlign { Absent, Malformed,
Parsed { blocks, end } }, computed only for same-contig mapped pairs (the
only reads any consumer uses it for), and thread &MateAlign through the
three callers:
  - mate_alignment_end -> mate_end_or_tlen(record, mate_start, &MateAlign):
    reads the precomputed end, same TLEN fallback for Absent/Malformed.
  - resolve_mate_span borrows the parsed blocks via Cow (owned only for the
    rare MC-absent TLEN estimate), so MateSpan no longer allocates on the
    exact path.

Absent conflates 'not a same-contig pair' with 'no MC tag', which is
sound: every consumer that could see the former already returns before
using the value (different-contig strand reads have proper=false, so the
mate end is unused). Malformed maps to the same skip the old  produced.

Output byte-identical on the 80M-read ENCODE BAM. Instructions retired
drop ~5.3% (~573.9B -> ~543.5B) and wall time ~29.2s -> ~27.3s this run;
cumulative with the two prior commits ~9.5% fewer instructions vs the
start of this pass.

rec_alignment_blocks, mate_alignment_blocks, and rec_introns each returned
a freshly heap-allocated Vec per read; on the 80M-read BAM that is on the
order of 150M+ small allocations, and rpmalloc alloc/dealloc was ~8% of
the compute thread's self time in the profile.

Return a SmallVec<[(u32,u32); 4]> (aliased CigarBlocks) instead: a read's
M/=/X blocks (unspliced 1, one junction 2, ...) and its introns fit inline
in the common case, so the parse never touches the heap; only reads
spanning more than four blocks spill. Every consumer already took a slice
(&[(u32,u32)]) or borrows via MateSpan's Cow, so nothing downstream
changed. The MateAlign::Parsed blocks field carries the same SmallVec.

Output byte-identical on the 80M-read ENCODE BAM. Instructions retired
drop ~3.75% (~544.0B -> ~523.6B) and wall time ~27.6s -> ~26.6s;
cumulative over this pass ~12.8% fewer instructions vs the start.

accept() allocated a fresh Vec<u32> for the gene loci overlapping each
read. Most reads fall in zero or one locus, so this was ~one small heap
allocation per read for a 1-2 element vector.

Return the sweep results into a SmallVec<[u32; 4]> (aliased LociVec)
instead: zero/one/few-locus reads stay inline and never touch the heap;
only reads over more than four overlapping genes spill. overlapping_loci
still clear()s and extend()s it, and every consumer already took &[u32],
so nothing else changed.

Output byte-identical on the 80M-read ENCODE BAM. Instructions retired
drop ~1.4% (~523.6B -> ~516.4B); cumulative over this pass ~14% fewer
instructions vs the start.

parse_mate_align went through mate_cigar, which did record.aux_tag(MC)?
.as_str().map(|s| s.to_vec()): aux_tag already hands back an owned
AuxValue (the aux store is cloned out), and .to_vec() then allocated and
memcpy'd a *second* copy of the same bytes for every paired read carrying
an MC tag.

Inline mate_cigar into parse_mate_align and move the Vec out of the owned
AuxValue::String via a let-else, so there's exactly one allocation (the
one aux_tag already made) and no extra memcpy. A missing or non-string MC
tag maps to Absent exactly as the old Option path did.

Output byte-identical on the 80M-read ENCODE BAM. Instructions retired
drop ~1.7% (~516.6B -> ~507.9B); cumulative over this pass ~15.4% fewer
instructions vs the start.

On the multi-locus classification path, classify_blocks unioned the
overlapping genes' exon and coding intervals by concatenating them and
calling merge_intervals, which sort_unstable()s the whole set. But each
locus's exon_union / coding_union is already sorted and disjoint by
construction, so the sort is redundant work — and quicksort was ~3% of
the compute thread's self time in the profile, almost all of it from this
per-read path.

Add union_of_sorted_disjoint (a pairwise merge of pre-sorted sets) and
use it for the coding and exon unions; the per-locus spans stay on
merge_intervals since there are only a few and they aren't pre-sorted as a
set. Factor the coalescing step both share into push_coalesced. The union
is the same set regardless of merge order, so output is unchanged; added
unit tests pin the k-way result to merge_intervals' and cover overlap,
adjacency, dedup, and empty/single-set inputs.

Output byte-identical on the 80M-read ENCODE BAM. Instructions retired
drop ~3.5% (~507.8B -> ~490.1B) and wall time ~26.2s -> ~24.8s;
cumulative over this pass ~18% fewer instructions vs the start.

Code review found that bases_overlapping_exons's merge-walk (binary-search
+ early-exit) silently assumed a transcript's exons are pairwise disjoint —
a precondition the parsers never enforce (they only sort). A malformed or
duplicate exon record for one transcript (from hand-edited or merged
annotation files) breaks the partition_point's monotonic-end precondition,
which could undercount to zero or, in debug builds, panic on u32 underflow.

Make the invariant real: normalize_exons() sorts each transcript's exons
and coalesces any that overlap, so the exon list handed to every downstream
interval walk is sorted and pairwise disjoint. On real annotations (whose
per-transcript exons never overlap) this only sorts — a no-op that keeps
output byte-identical on the 80M-read ENCODE BAM. It does not change the
locus exon_union at all (merge_intervals already coalesced the same
intervals), so gene-level coverage is untouched; the only effect is that a
malformed overlapping-exon transcript now counts each exonic base once
(sane deduplicated per-transcript coverage) instead of double-counting.

Keeps the merge-walk's full speed (the ~2.3% it's worth); normalize is a
one-time load cost (~0.3%).

…d-MC coverage

Follow-ups from the perf-pass review, no behavior change:
- Delete a stale doc line left over from the removed mate_cigar fn, which
  had glued itself onto parse_mate_align's doc block.
- Move the smallvec import into the third-party group (was appended after
  the crate:: group).
- Rename mate_alignment_end_includes_trailing_skip -> _blocks_ (it exercises
  mate_alignment_blocks; the old function name no longer exists).
- Point the bases_overlapping_exons doc at the newly-enforced disjoint-exon
  invariant.
- Add insert_size_skips_pair_with_malformed_mc: pins that a present-but-
  unparseable MC tag (MateAlign::Malformed) skips the pair rather than
  falling through to the TLEN estimate.